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ABSTRACT 

Fourier spectral method can achieve exponential accuracy both on the approximation level 
and for solving partial differential equations if the solutions are analytic. For a linear partial 
differential equation with a discontinuous solution, Fourier spectral method produces poor point- 
wise accuracy without post-processing, but still maintains exponential accuracy for all moments 
against analytic functions. In this note we assess the accuracy of Fourier spectral method applied 
to nonlinear conservation laws through a numerical case study. We find that the moments with 
respect to analytic functions are no longer very accurate. However the numerical solution does 
contain accurate information which can be extracted by a post-processing based on Gegenbauer 
polynomials. 
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1 Introduction 


In this note we are concerned with the accuracy of Fourier spectral method for the solution of a 
nonlinear conservation law 


d t u + d x f(u) 
u(x , 0) 


0 , - 1<*<1 

u°(x) 


( 1 . 1 ) 


where the initial condition u°(x) is 2-periodic. As is well known, solutions to (1.1) typically con- 
tain discontinuities even if the initial condition u°(ar) is analytic (In this paper, for simplicity of 
presentations, we will use analytic functions to represent general smooth functions; similar results 
can also be obtained for C k or C°° functions). The purpose of this note is to assess accuracy under 
such a situation through a numerical study. 

We start by recalling the Fourier approximation operator Sat to an L 2 function u(x): 

N 

S N u(x)= £ u k e ik ™ (1.2) 

k=-N 

where the Fourier coefficients Uk are defined by 

u k = ^j\(x)e- ik ™dx (1.3) 


for Fourier Galerkin, and by 

= 2NTT £ “ (i ’ k 


N 


tk7TXi 


2i 


j=-N 


Xj 2 N + 1 


(1.4) 


for Fourier collocation. We will also use the notation uj^(x) = Sjsf ti(z). To solve the partial 
differential equation (1.1), the standard Fourier spectral algorithm is 


SN(dtVN + dxf{vN )) = 0, -1 < x < 1 

vn(x,0) = Snu°(x) (1.5) 

where viy(x,t) = v k (t)e tkirx is supposed to approximate the exact solution u(x,t) of (1.1), 

and 5jv is either the Galerkin or the collocation Fourier approximation operator defined by (1.2)- 
(1.3) or by (1.2)-(1.4). 
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The approximation error 


u(x) - Snu(x) 


( 1 . 6 ) 


is well known to be exponentially small (i.e., it is 0(r N ) for 0 < r < 1) if u(x ) is analytic. However, 
if u(z) is only piecewise analytic and discontinuous, the approximation error (1.6) is 0(1) near the 
discontinuity and 0(jj) elsewhere. This is known as the Gibbs phenomenon (see, e.g., [4] and [3]). 
Fortunately, even if the accuracy is poor in the point-wise sense, it is still excellent for the moments 
against any analytic functions. For any L 2 function n(x) and any analytic function w{x), we have 
[5]: 


J (n(x) — uj^(x))w(x)dx 


< Cr 


.N 


(1.7) 


for some constant C and 0 < r < 1. This property is the basis of all the “reconstruction” or “post- 
processing” techniques. These techniques try to recover exponential or at least high order accuracy 
for point values based on the Fourier approximation S^u(x) of a piecewise analytic function. In 
other words, one tries to obtain a small post-processed approximation error 


u(x) - PjvSjvu(x) (1.8) 

where P /v is some post-processing operator. Examples of P^ include various high frequency filters 
[14], [11], [16], [2], which are of the form 

N / k\ 

P N S N u(x)= ^2 a ( — ) u k e' kirx (1.9) 

k=-N ' 

with Sjvu(x) given by (1.2). The function cr(f) in (1.9) is even (or satisfies cr(— £) = <r(f) if it 
is complex valued as in [2]), smooth (the accuracy of the filter depends upon its smoothness), 
supported in (-1,1) and satisfies <t(0) = 1 and a^(0) = 0 for 1 < k < K (with accuracy of the 
filter again depending on K ). These filters can recover high order or even exponential accuracy 
in the smooth regions away from the discontinuities (the filter in [2] can also recover high order 
accuracy up to the discontinuity from one side). A more recent example of Pjy is the Gegenbauer 
polynomial based procedure discussed in [6], [7], [8], [9] and [10], which can give uniform exponential 
accuracy for all x right up to the discontinuity for piecewise analytic functions. In this sense spectral 
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Fourier approximation is also exponentially accurate for piecewise analytic functions — one only has 
to extract the hidden information from the poor approximation Sn(x) using the post-processor 

Pn . 

When spectral method is used to solve the PDE (1.1), we can consider the following different 
types of errors. The strongest is the point-wise error from the exact solution u(z,f): 


u(x,t) - v N (x,t), (1.10) 

which cannot be small even for t = 0 due to the Gibbs phenomenon. A more reasonable error is 
the point-wise error of the numerical solution Vjv(x 7 t) from the Fourier approximation of the exact 
solution ujv(x, t ): 

u N (x,t) - v N (x,t). (1.11) 

If this error is exponentially small, we claim the spectral method for (1.1) is exponentially accurate 
because of the post-processor (1.8) for the exact solution u(x,2). An even weaker error is defined 
by the error in the first few Fourier coefficients, i.e. 


Uk(t) - v k (t) 


( 1 . 12 ) 


for the first few fc, where Uk(t) are the Fourier coefficients of the exact solution u(a:,<) of (1.1). 
This is actually an example of the more general definition of error in moments with respect to any 
analytic function w(x): 

j (un(x) - vn(x))w(x)cLx (1.13) 

In fact, as long as this error in moments is exponentially small, we claim that the spectral method 
is exponentially accurate in solving (1.1) by using property (1.7) for the exact solution u(x,t) and 
the post-processing (1.8). 

If the PDE (1.1) is linear (i.e. f(u ) = a(z,t)u), it is proven in [5], [1] that spectral Fourier 
method is exponentially accurate in the sense that (1.13) is exponentially small. A post-processing 
(1.8) applied to vjy(x, £) would then yield an exponentially accurate point-wise approximation to the 
exact solution u(x y t). However, if (1.1) is nonlinear, it is still a theoretically open problem whether 
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spectral Fourier method, equipped with either high frequency filtering or vanishing viscosity [15], 

[12] , is exponentially (or high order) accurate in the sense of (1.13). Computational evidence in 

[13] seems to suggest that, even in this nonlinear case, highly accurate information is still implicitly 

contained in the numerical solution and can be extracted (at least away from the discontinuity) 

by a post-processing using high frequency filtering. In the next section we will perform a detailed 

2 

numerical case study about this accuracy issue for Burgers’ equation ( f{u ) = We use a 

high frequency solution filter to stablize the algorithm, and post process the numerical result 
using the procedure based Gegenbauer polynomials [6], [10]. We observe that the spectral Fourier 
method is not very accurate in the sense of moments against analytic functions (1.13). However, 
numerical evidence does indicate the possibility of very high accuracy under some weaker definition 
of accuracy, perhaps some average of Fourier coefficients, since the post-processed result Pnvn{x, t) 
is much more accurate than the Fourier coefficients themselves, and accurate Fourier coefficients 
can be “reconstructed” from this post-processed solution Pjvvn(x, t). 

2 A Numerical Case Study about Accuracy 

In the numerical solution reported in this section, time discretization is by a third order Runge- 
Kutta method, with a time step At sufficiently small such that the spatial error is dominant 
in all cases. We compute the exact solutions of the PDE by Newton iterations on the implicit 
characteristic equations, and compute the Fourier coefficients of a function (if not analytically 
given) by using a sufficiently accurate numerical quadrature. 

We first solve a linear equation 

= 0, -1 < x < 1 

= X (2.1) 

with periodic boundary conditions, up to t=l, using the Fourier Galerkin method: 

Sn (d t VN + - — — - 8 x vn) = 0 , 

\ 5 — 4cos(7t:e) / 


d t u + 


(-);r Xi 


5 — 4 cos(7rx) 

u(x, 0) 
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( 2 . 2 ) 


N (-U k i 

v N (x, o) = s NX = y, ~ih eikrx 

k = -N 

0 


Standard Galerkin method is stable for this linear problem but produces poor point value accuracy 
(Figure 1, left). However, the accuracy in the first few Fourier coefficients, as representatives of 
moments against analytic functions, are computed more accurately (Figure 1, right). 



Fig. 1: Errors in log scale, linear PDE (2.1). Fourier Galerkin using 2 N + 1 modes, for N — 10; 
N = 20; N = 40 and N = 80. Left: point-wise errors; Right: errors in the first 10 Fourier 
coefficients. 

In order to compare with the nonlinear case reported later, we solve the same linear equation 
(2.1) using the filtered Fourier method, i.e., after each Runge-Kutta time step, the numerical 
solution is filtered by (1.9) with the exponential filter: 

<r(0 = (2.3) 

where r is increasing with N and is related to the order of the filter, and c* is chosen such that e a 
equals machine zero ( 10~ 16 for double precision). The exponential filter (2.3) has the advantage 
of simplicity, and usually it works equally well as more complicated filters [16]. For this linear 
problem, as well as for the nonlinear Burgers’ equation later, we will use the Fourier method with 
the following choice of filter orders: r — 4 for N — 10; r = 6 for N = 20; r = 8 for N = 40 and 
r = 12 for N = 80. The result is shown in Figure 2. Comparing with Figure 1, we can see better 
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point value accuracy in the smooth region because of the filters, and similar (good) accuracy for 
the first few Fourier coefficients. 



Fig. 2: Errors in log scale, linear PDE (2.1). Fourier Galerkin using 2 N + 1 modes with exponential 
solution filters of order r. r = 4 for TV = 10; r = 6 for N = 20; r = 8 for N = 40 and r = 12 for 
N = 80. Left: point-wise errors; Right: errors in the first 10 Fourier coefficients. 

The computational result for the linear equation is not surprising since it just shows the proven 
fact [5], [1] that Fourier coefficients, as representatives of moments against analytic functions, are 
computed with exponential accuracy by the spectral Fourier method, and filtering will recover 
exponential point value accuracy in smooth regions away from the discontinuity. It should be 
noticed that, for the same N y the accuracy for the first few Fourier coefficients is at the same level 
or better than the best point value accuracy in the smooth region after filtering. This is again not 
surprising since point value accuracy is obtained from the Fourier coefficients through filtering. 

We now consider the nonlinear problem we are really interested in: we solve the nonlinear 
Burgers’ equation 

= 0, — 1 < x < 1 

= 0.3 4- 0.7 sin(7rz). (2.4) 

The solution develops a shock at t — and we compute the solution up to t = 1. The initial 
condition is chosen such that the shock is moving with time. For this nonlinear PDE, the standard 


d t u + d x ( — 


u(x , 0) 
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Galerkin method cannot converge to the entropy solution [15]. One would need to add dissipations 
either by the high frequency solution filtering (1.9) or by the spectral vanishing viscosity [15], [12], 
[13]. Numerical results for the Burgers’ equation with the vanishing viscosity method can be found 
in, e.g., [13]. Here we will only report the results obtained by solution filtering, using the same r 
as in the previous linear case (2.1). We have also computed with the vanishing viscosity methods 
and have obtained similar results. 

In Figure 3 we plot the point-wise error u(x,t) — vw(x,t) (left), and the error for the first 10 
Fourier coefficients (right). While the pattern of the point-wise errors are similar to the Unear case 
in Figure 2, the errors for the Fourier coefficients are clearly much worse in comparison. As a matter 
of fact, for the same TV, the errors for the first few Fourier coefficients are a few magnitudes larger 
than the smallest point value error in the smooth region. This is clearly different from what we 
observe in the Unear case in Figure 2, and suggests that the first few Fourier coefficients, again as 
representatives of moments against analytical functions, are no longer computed with exponential or 
high order accuracy. It is sort of puzzhng that each difference in the Fourier coefficients 
is relatively large (Figure 3 right), but the point-wise error u(x, t) — visi(x, £), which is just an average 
of Uk(t) — Vk(t) (against an 0(1) function e tklTX ), is much smaller in the smooth region (Figure 3 
left). Clearly some canceUation is present. 



Fig. 3: Errors in log scale, Burgers equation (2.4). Fourier Galerkin using 2 N + 1 modes with 
exponential solution filters of order r. r = 4 for N = 10; r = 6 for N = 20; r = 8 for N ~ 40 and 
r — 12 for N = 80. Left: point-wise errors; Right: errors in the first 10 Fourier coefficients. 
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Next, we apply the Gegenbauer post-processor [6] to v^(x^t). This procedure can be roughly 
described as follows: given the Fourier partial sum uj^(x) of an analytic but not periodic function 
w(ar), one first finds the approximations to the first m Gegenbauer expansion coefficients of the 
function u(x). Here Gegenbauer polynomials are orthogonal polynomials in [—1, 1] under the weight 
function (1 — z 2 ) A ~ 2 . One then uses this Gegenbauer series with the first m terms to approximate 
u(x) everywhere in [—1,1]. To use this procedure, one must know the location of the discontinuity 
(however, the procedures in [8] allows one to handle the case where the location is not known 
exactly), and to choose the parameters A and m. It is proven in [6] that when m and A are both 
chosen proportional to (but less than) TV, the reconstructed point values are exponentially accurate 
everywhere inside [—1,1]. Thus Gibbs phenomenon is completely removed. The details can be 
found in [6], [7], [8], [9] and [10]. 

We would like to point out that there is no theoretical justification in doing this post-processing 
for the current nonlinear case, since the post-processing procedure assumes that the Fourier coef- 
ficients are accurate, which is not true any more. However, the post-processed result is surpris- 
ingly good (Figure 4). We can observe good accuracy everywhere including at the discontinuity 
x = ±1 + 0.3 as in the approximation test case discussed [6]. The reconstructed Fourier coef- 
ficients, namely the Fourier coefficients of Pj^v^(x y t), are much more accurate than before the 
post-processing (compare Figure 4 right with Figure 3 right). 



Fig. 4: Errors in log scale, Burgers equation (2.4). Fourier Galerkin using 2N + 1 modes with 
exponential solution filters of order r. r = 4 for N — 10; r = 6 for N = 20; r = 8 for N = 40 
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and r = 12 for N = 80. Gegenbauer post-processed, with parameters A = 2 ,m = 1 for N = 10; 
A = 3,m = 3 for N = 20; A = 12, m = 7 for JV = 40 and A = 62, m = 15, for N = 80. Left: 
point-wise errors; Right: errors in the first 10 Fourier coefficients. 

This suggests that, even if or its Fourier coefficients are not very accurate, it 

contains accurate information which is extracted in this case by the Gegenbauer polynomial based 
post-processor ffy* This numerical evidence suggests that in the nonlinear PDE case, Fourier 
coefficients £*(/), just like point-wise values in the linear (or nonlinear) PDE case, are no longer good 
indicators of accuracy. They themselves are not very accurate, but they implicitly contain accurate 
information which can be extracted by adequate post-processors Pjy. This accurate information 
might be contained in some averages of the Fourier coefficients (since the post-processing procedure 
based on Gegenbauer polynomials [6] uses certain averages of Fourier coefficients rather than the 
coefficients themselves). 

We finally make two remarks: 

Remark 2.1. In the Gegenbauer reconstruction procedure above we have used the exact shock 
location. The procedure in [8] allows us to use an approximate shock location, determined from 
the Fourier coefficients themselves (e.g., [2]). Similarly good results can be obtained when the 
reconstruction is performed in a slightly smaller sub-interval inside which the solution is guaranteed 
to be analytic. For example, we use the shock detector in [2], which in this case detects the shock 
location to within 0.0000025 for all the N values used, and a reconstruction inside the sub-interval 
[—0.999997,0.999997], which is just slightly smaller than [—1,1] (when numerically detected shock 
is shifted to x = — 1) and guarantees that the true shock is outside this region. The result is shown 
in Figure 5 (left). It is clearly as good as the case where the exact shock location is used (compare 
with Figure 4 left). 

Remark 2.2. If we use collocation (1.4) instead of Galerkin, (for the reconstruction procedure, 
see [10]), the result is almost identically good: Compare Figure 5 (right) with Figure 4 (left). 
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Fig. 5: Point-wise Errors in log scale, Burgers equation (2.4). Fourier method using 2N + 1 modes 
with exponential solution filters of order r. r — 4 for N = 10; r = 6 for N = 20; r ~ 8 for N = 40 
and r = 12 for N = 80. Left: Galerkin, Gegenbauer post-processed with a numerically determined 
shock location using the techniques in [2], which for this problem produce shock locations to within 
0.0000025 for all the N used. The reconstruction sub-interval is [—0.999997,0.999997] when the 
numerical shock is shifted to a; = -1. Parameters: A = 2,m = 1 for N = 10; A = 3,m = 3 
for N = 20; A = 26, m = 9 for N = 40 and A = 52, m = 17, for N — 80. Right: collocation. 
Gegenbauer post-processed, with parameters A = 2, m = 1 for N = 10; A = 3, m — 3 for TV = 20; 
A = 26, m = 9 for TV = 40 and A = 60, m = 15, for N = 80. 

3 Concluding Remarks 

Through a careful numerical case study for the Burgers’ equation, we have found that the Fourier 
spectral method, equipped with spectrally small dissipations in the form of high frequency filters 
or vanishing viscosities, are not accurate in the first few Fourier coefficients, or in moments against 
smooth functions. However, accurate information is indeed contained in the numerical solution, 
and can be extracted by using the Gegenbauer polynomial based post-processor in [6]-[10], 
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